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We discuss energy densities in the strong-interaction limit of density functional theory, deriving 
an exact expression within the definition (gauge) of the electrostatic potential of the exchange- 
correlation hole. Exact results for small atoms and small model quantum dots are compared with 
available approximations defined in the same gauge. The idea of a local interpolation along the 
adiabatic connection is discussed, comparing the energy densities of the Kohn-Sham, the physical, 
and the strong-interacting systems. We also use our results to analyze the local version of the Lieb- 
Oxford bound, widely used in the construction of approximate exchange-correlation functionals. 



I. INTRODUCTION 

Increasing the accuracy of the approximations to the 
exchange-correlation energy functional E xc [p] of Kohn- 
Sham (KS) density functional theory (DFT) is of cru- 
cial importance for research areas ranging from theoreti- 
cal chemistry and biochemistry to solid-state and surface 
physics (for a recent review, see, e.g., Ref. [T]). 

A piece of exact information on E xc [p] is provided by 
the strong-interaction limit of DFT, in which the cou- 
pling constant of the electron-electron interaction be- 
comes infinitely large while the one-electron density p(r) 
does not changeP^ This defines a fictitious system with 
the same density of the physical one and maximum possi- 
ble correlation between the relative electronic positions, 
useful to describe situations in which restricted Kohn- 
Sham DFT encounters problems, such as low-density 
man y-pa rticle scenarios and the breaking of the chemical 
bond! 6 -' 7 -' The exact mathematical structure of this limit 
has been uncovered only recentlyP^ anc j exact calcula- 
tions (at least f or simple systems) have started to become 
available! 6 * 7 * 11 * 

The aim of this paper is to make a step forward in 
the inclusion of this new piece of exact information into 
approximations to the exchange-correlation energy func- 
tional E xc [p\. Previous attempts in this directiorPHSl fo- 
cused on global (i.e., integrated over all space) quantities, 
introducing size-consistency errors. The exact solution of 
the strong-interaction limit, now available, makes accessi- 
ble not only global, but also local quantities, from which it 
is easier to construct size-consistent approximations^^ 
(for a critical review on size consistency of approximate 
energy density functionals see, e.g., Ref. [15] and, espe- 
cially, Ref.QH). 

Local quantities are in general not uniquely defined 
(for a complete discussion see, e.g., Ref. IT?]) . Here we 
focus on the conventional, physical transparent, defini- 
tion in terms of the electrostatic energy of the exchange- 
correlation hole. We derive the exact expression for this 
energy density in the strong-interaction limit, and we 
evaluate it for small atoms and small model quantum 
dots, making comparisons with available approximations 
within the same definition (same gauge). We then dis- 



cuss the idea of a local interpolation along the adiabatic 
connection by comparing energy densities in the physical 
case, the weak- and the strong-interaction limit. 



As a byproduct, our results allow to analyze the 
local version of the Lieb-Oxford bound, a condi- 
tion widely used to construct approximate exchange- 
correlation functionals. As well known, the Lieb-Oxford 
bound is an exact conditiorP^H22l on th e global E xc [pj. 
Many non-empirical approximate functionals, however, 
use its local version, which is a sufficient but not nec- 
essary condition to ensure the global bound (see, e.g., 
Refs. 1211221 and 123 ) . Our analysis strongly suggests that 
the local version of the Lieb-Oxford bound is formulated 
in the gauge of the electrostatic potential of the exchange- 
correlation hole. This, in turn, implies that the local 
bound is certainly violated at least in the tail region of a 
molecular or atomic density, and in the bond region of a 
stretched molecule. 



The paper is organized as follows. In the next Sec. |TT] 
we review the DFT adiabatic connection as a tool to build 
approximate E xc [p], highlighting the role of the strong- 
interaction limit and discussing the size-consistency prob- 
lem of interpolations based on global quantities. In 
Sec. |III| we discuss energy densities in general and we in- 
troduce the gauge of the electrostatic energy associated 
to the exchange-correlation hole. The exact expression 
for this quantity in the strong-interaction limit is derived 
in Sec. |IV| where also approximations are discussed. In 
particular, we analyze the "point-charge plus continuum" 
(PC) model functional,^ showing that it is meant to be 
an approximation to the energy density within the same 
conventional definition considered here. Energy densi- 
ties along the adiabatic connection are discussed and an- 
alyzed in Sec. |Vj We then use our results to discuss the 
local version of the Lieb-Oxford bound in Sec. IVII The 
last Sec. |VIl| is devoted to conclusions and perspectives. 
Finally, a simple illustration for the case of the uniform 
electron gas is reported in the Appendix. 
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II. INTERPOLATION ALONG THE 
ADIABATIC CONNECTION 

Within the framework of the adiabatic connectiorP^H 
the exchange-correlation energy can be expressed by the 
coupling constant integration 

E xc [p}= f d\(y x [p}\V ee \y x [p})-U[p] (1) 
Jo 

r-l 

dXW x [p], (2) 

'o 

with ^ \[p] being the ground state wave function of a 
fictitious system with scaled electron-electron interaction 



H = T + XV e , 



V a. 



(3) 



The external potential V^ xt is adjusted to keep the den- 
sity p\(r) out of ^>\[p] in agreement with the physical 
density, p\(r) = pi(r) = p(r). For the weak interaction 
limit (A = 0) we encounter the Kohn-SharrP^l reference 
system and the integrand Wo[p] becomes the exact ex- 
change energy E x [p\. For the strong-interaction limit 
(A — > oo) a reference system within the strictly corre- 
lated electrons concept (SCE^EMU can be defined. The 
asymptotic expansions of W\ [p] are 



Wx^o[p] = E x [p}-\ 

W X ^oo[p} = W 00 [p] 



2A£ c GL2 [p] 

, w L[p\ , 



+ 0(A 2 ) 
0(A- p ), 



(4) 
(5) 



where E^ h2 [p] is the correlation energy given by second 
order Gorling-Levy perturbation theory (GL2j2Sand p > 
5/4 P Exact expressions for the functionals W^oofp] and 
[p] are given, respectively, in Refs. [8| and [9l 

Expression ([I]) for the exchange-correlation energy is 
exact as lo ng as the exact dependence of the integrand on 
A is known J22G2D As this is obviously not the case, Eq. 
still enables approximate exchange-correlation energies 
by modeling W\ [p] along the adiabatic connection. 

Attempts towards approximate Wa[p] have been un- 
dertaken by BeckeP^ introducing the half and half func- 
tional. A model is defined assuming a linear dependence 
of W\ [p] on A and setting Wq [p] equal to exact exchange 
and to LSDA exchange-correlation. This results 

in a functional with 50% exact exchange and 50% LSDA 
exchange-correlation. Further adjustment of the portion 
of exact exchange by semi-empirical arguments gives rise 
to hybrid functionals like B3LYPp3MlThe adiabatic con- 
nection may also be used for the construction of non- 
empirical hybrids as in Ref. 36. Here a model for WaIp] 
is defined consisting of two intersected straight lines fixed 
by exact exchange, GGA exchange and GGA exchange- 
correlation. ErnzerhoP-^ introduced a curved model by 
proposing a Pade interpolation for the integrand using as 
input exact exchange and E^* L2 [p] in the weak interaction 
limit and GGA exchange-correlation for A = 1. 

The mentioned models for the integrand (i.e. all except 
B3LYP) have in common that for the weak interaction 



limit exact exchange is employed and for the physical sit- 
uation with A = 1 information from approximate DFT 
(DFA) is used. The argument for the recourse to exact 
exchange is that DFA exchange works well only if com- 
bined with DFA correlation. This is due to error cancel- 
lation. Consequently, DFA exchange-correlation can be 
used for the physical case where exchange and correla- 
tion are employed together. As error cancellation in DFA 
exchange-correlation might not be satisfactory, a contin- 
uation of the ansatz of Ernzerhol^ is possible by taking 
DFA exchange-correlation at some intermediate A instead 
of A = 1. This would allow to balance the exchange error 
with the correlation error. Along this line, Mori-Sanchez, 
Cohen and Yang 3 ^ constructed their MCY1 functional: 
a Pade interpolation is undertaken with exact exchange 
and meta-GGA exchange input in the weak interaction 
limit and meta-GGA exchange-correlation for an inter- 
mediate A (chosen semi-empirically) . 

The discussed models clear ly ou tperform the stand 
alone DFAs they are based on! 31 | 36| E^ Nonetheless, em- 
ployment of DFA quantities in their construction can lead 
to serious misbehavior in the curvature of the integrand 
as demonstrated by Peach, Teale and Tozer^by compar- 
ison of the MCY1 approximation with accurate quanti- 
ties along the adiabatic connection (see, e.g., Fig. 3 in 
Ref. l39|) . In the same paper the authors show that accu- 
rate exchange-correlation can be recovered by establish- 
ing an interpolation with accurate full-CI ingredients. 

A model that avoids unfavorable DFA bias is the inter- 
action strength interpolation (ISI) PSED Here information 
from the weak interaction limit is employed, namely ex- 
act exchange and GL2, together with information from 
the strong interacting limit through the point-charge 
plus continuum (PC) model, which provides approxi- 
mate expressions for W^p] and W^[p\. The A depen- 
dence of Wa[p] is then modeled by an interpolation be- 
tween the two limits. Nowadays, the functionals W^oofp] 
and W^p] can be accur ately computed within the SCE 
many-electron formalismP^l Refs. [5] and [§] compare the 
PC solutions with the exact SCE values for small atoms: 
while W^. [p] is a very reasonable approximation to its 
exact counterpart P the original W^ c [p] turned out to be 
much less accurateP The exact results could be used to 
propose a revised PC approximation W^ vFC [p] having 
accuracy similar to the one of [p]. Further com- 
parison is undertaken in section [IV B| of this paper for a 
more refined quantity, the local energy density as it will 
be defined in the next sections. 

Although unplugged from any DFA bias (if we use ex- 
act input quantities), an unpleasant feature of the ISI is 
the violation of size consistency. This is due to the non- 
linear way the (size-consistent) ingredients Wo[p], ^oMi 
VFoofp] and W^lp] enter the interpolation. For example, 
the revised ISI (which behaves better in the A — > oo limit 
than the original ISI) reads^ 



wr Isl [ P ] 



a[p]X + 



b[p}\ 



y/l + c\p]X + d[p] 



(6) 



3 



where a, b, c and d are non linear functions of Wo[/j], 
Wq[/o], Woo[p] and W^Jp], determined by imposing the 
asymptotic expansions of Eqs. Q and 



In the definition of Eq. (12 1, u>A[p](r) is the electrostatic 



a[p] 
b[ P ] 

c[ P ] 



d[p] = -1- 



(E x [p\ - WooM) 2 
{E x \p]-W 00 [p]) i 

ze^pW^p? 



{E x \p]-WM)*' 



(7) 
(8) 

(9) 
(10) 



Notice that the lack of size consistency is shared by all 
functionals in which the exact exchange energy (or any 
global energy) enters in a nonlinear way. Thus also, for 
example, MCY1. 

As a final remark on the revISI functional, we can add 
that if one makes the approximation E^ L2 w E^f P2 , it 
can be viewed as a double hybrid functional (see, e.g., 
Refs. l40"H43|) . With respect to available double hybrids, 
the revISI lacks size consistency, but it has the advantage 
of being able to deal with the small-gap systems problem- 
atic for perturbation theory. The practical impact of the 
lack of size consistency of the revISI functional still needs 
to be tested. 



III. ENERGY DENSITIES: DEFINITIONS 

A possible way to recover size consistency in the ISI 
framework is to use a local integrand in Eq. ([!]) : 



E xc [p] 



dX / drp(r) w x [p](r). 



(11) 



The idea is then to build a local model, u^ SI [p](r), by in- 
terpolating between the A — > and the A — > oo limits. As 
said, the energy density w\[p)(r) is not uniquely defined, 
so that an important requirement here is that the input 
local quantities in the weak- and in the strong-interaction 
limits are defined in the same way (same gauge). 

One of the most widely used definitions of the energy 
density in DFT is in terms of the exchange-correlation 
hole (see, e.g., Refs. QUEH and HD h* c (r,r'), 



where 



. If /4(r,r') , 



^ c (r,r') = ^%^-p(r'), 



Pir) 



(12) 



(13) 



and the pair-density PHj, ^) is obtained from the wave 
function * A [p] of Eqs. Q-Q: 

P 2 A (r,r') = N(N — 1) x 
^2 J l*A(roi, r'cr 2 , r 3 cr 3 , r N a N )\ 2 dr 3 ...dr N . (14) 



potential of the exchange-correlation hole (a negative 
charge distribution normalized to —1), around a refer- 
ence electron in r. This quantity at the physical coupling 
strength A = 1 (plus the Hartree potential) has been also 
called Vcond(r) in the literature (see, e.g., Ref. H1J). 

The energy density w\ [p] (r) in the A — > oo limit in the 
gauge of Eq. ( 12 ) is the central quantity of this paper: 



we will derive in the next section an exact expression 
using the strictly-correlated electron concept, and we will 
evaluate it for small atoms and quantum dots. Notice 
that the relevance of ^^[^(r) for constructing a new 
generation of approximate E xc [p] has also been pointed 
out very recently by BeckeP^ 



IV. ENERGY DENSITIES IN THE STRONG 
INTERACTION LIMIT 

A. Exact 

When A 4 oo the wave function ^\[p] tends to the 

strictly-correlated electron state, ^a^ooIp] — > *sce[p], 
witrPS 



|*SCE(ri,...,rjv)| 2 = 



1 

v 



ds^(ri-Ws)) 



x6(r 2 - f P(2) (s))...(5(r A r - fp (Ar) (s)), 



(15) 



where fi, ..,fjv are "co-motion functions", with fi(r) = r, 
and V denotes a permutation of {1,...A}. This means 
that the N points ri,...,rjy in 3D space found upon si- 
multaneous measurement of the N electronic positions in 
the SCE state always obey the N — 1 relations 



r, = f((ri) (i = 2,...,N). 



(16) 



In other words, the position of one electron determines all 
the relative N—l electronic positions (limit of strict cor- 
relation). All the N—l co-motion functions f,(s) satisfy 
the differential equation 



p(fc(r))dfc(r) = P(r)dr, 



(17) 



which, together with the group propertie d 1 1 of the (r) , 



ensure that the SCE wavefunction of Eq. (15) yields the 



given density p(r). Equation (17) has also a simple phys- 



ical interpretation: since the position of one electron de- 
termines the position of all the others, the probability 
of finding one electron in the volume element dr about r 
must be the same of finding the i th electron in the volume 
element dfj(r) about f,(r). 

Inserting Eq. (15) into Eq. (14) we obtain for the 
pair density P 2 A_S 
interaction limit 



i"i,r 2 ) = P 2 3( ^r 1 ,r 2 ) in the strong- 



P| CE (r l5 r 2 ) = |; / 



A ^( S ) XI 

— <5(ri 



f l (s))<5(r 2 -f J (s)), 



(18) 
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which also has a transparent physical meaning: two elec- 
trons can only be found at strictly correlated relative 
positions. 



We first compute 
P 2 SCE (r,r') 



ds 



N 



(19) 



where, in the right-hand-side, we have already inte- 
grated over the variable r'. From the properties of the 
Dirac delta distribution and of the co-motion functions, 
Eq. Jl9|) becomes 



ff CE (r, r') 



dr' 



1 



JV 

N ^ 



p(fr 1 (r)) \detd a f^(r)\ 



MS 



P(r) 

AT 



JV 



Slr-f^f-W)!' (20) 



where |det 9 a /j ^ (r)| (with a,j9 = x, y, z) is the determi- 
nant of the Jacobian of the transformation r — > f,~ (r), 
and we have used the fact that all the f»(r) (and their 
inverses, which, by virtue of the group properties of the 
co-motion functions are also co-motion functions for the 
same configuration^!) satisfy Eq. (17). Now we can 
use once more the group properties of the co-motion 
functions to recognize that for all i ^ j the function 
fj (f~ 1 (r)) must be another co-motion function with the 
exclusion of fi(r) = r (the identity can arise only if i = j). 
The double sum in the last term of Eq. ( 20 ) is then ex- 



actly equal to N times a single sum over all the co-motion 
functions ffc(r) with k > 2, so that 



^f CE (ry; 



dr' 



N 1 



k=2 



ffc(r) 



(21) 



Inserting Eq. (21 1 into Eqs. ( 12 )-( 13 ) we finally obtain 
1 



Woo(r) 



N 



1 



1 



k=2 



ffc(r) 



(22) 



where Vjj (r) is the Hartree potential. Notice that in pre- 
vious work the exact Woo [p] was given as^D 



U[p], (23) 



suggesting a corresponding energy density 



Woo(r) 



1 N 

-Y 



i=l 



2^|f i (r)-f,-(r)| 2 



/ 

(24) 



Equations ( 24 1 and ( 22 ) yield the same Woo [p] when in- 



tegrated with the density p(r), but are locally different. 
They show a general feature of the co-motion functions: 
any given energy density u>2o(r) can be always trans- 
formed into a different energy density (r) defined as 



,(r) 



1 



N 

-Y>. 

i=l 



(25) 



When multiplied by the density p(r), w^(r) and - u;^ (r) 
integrate to the same quantity, because all the co-motion 
functions (and their inverses, which are also co-motion 
functions) satisfy Eq. (17). Only Eq. (22) corresponds 
to the gauge of the exc 
Eqs. (pl-pih. 



lange-correlation hole defined by 



B. Approximations: the PC model 

The point-charge-plus-continuum (PC) modeP^I [ s a 
physically sound approximation to the A — > oo indi- 
rect electron-electron repulsion energy Woo[p]. The idea 
is to rewrite the indirect Coulomb interaction energy 
W\[p] as the electrostatic energy E es \^>\,p] of a sys- 
tem of N electrons in the state ^\[p] embedded in a 
smeared background of positive charge p+(r) = p(r)W 
In fact, this total electrostatic energy E es [$>\,p] is just 
the sum of the electron-electron repulsion energy, E ee — 

\\V ee \^ \) , the electron-background attraction energy, 
E eb = —2U[p], and the background-background repulsion 
energy E bb = U[p], thus yielding exactly E es [^\,p] = 
E ee + E Eb + E bb = W x [p}. 

This relation is valid for every A, but in the A — > oo 
limit, when ^>x —> ^sce, we expect that the electrons 
minimize E es [$>xi p] by occupying relative positions that 
divide the space into neutral cells with possibly zero 
(or weak) lowest-order electrostatic multipole momentsP 
The idea is then that for one of the SCE configura- 
tions {r, faO"), fzv(r)} we may approximate the indi- 
rect electron-electron repulsion by the sum of the electro- 
static energies of all the cells (i.e., we neglect the cell-cell 
interaction in view of their neutrality and low multipole 
moments): 

N 

e es (r, f 2 (r), f N (r)) a £ E ceQ ([p\; f^r)), (26) 
f=i 

where E cell ( [p] ; i\ ) is the electrostatic energy of the cell 
around an electron at position r<, equal to the sum of 
the attraction between the electron and the background 
contained in the cell and the background-background re- 
pulsion inside the cellP 

Notice that, for a given SCE configuration, the elec- 
trostatic energy e es (r, fa(r), fjv(r)) of Eq . (|26| is equal 
to Nwozfc), where Woo(r) is given in Eq. (|24[). The PC 
model is then trying to approximate Nw^\r) by con- 
structing the electrostatic energy E cc u([p\; r$) of a cell 



5 



around the electron at position r^. However, and this is 
a crucial step to understand the gauge of the PC model, 
once an approximation for _E cc n([p]; r^) has been built, 
the sum over the N electrons in the right-hand side of 
Eq. (26 1 is replaced by NE ce u ( [p] ; r) P In the original 
derivation of the PC modeP this step was seen as a fur- 
ther approximation. We now know, thanks to the exact 
SCE formulation, that this is is not an approximation, 
but an exact feature of the A — > oo limit, clarified in 
Eq. ( 25 ) . Because of this transformation, the local elec- 



trostatic energy that the PC model is trying to approxi- 
mate is then exactly the same as the one of the exchange- 
correlation hole of Eq. ( 22 ) . 



It is important to stress that the PC cell is not an 
approximation to the exchange-correlation hole in the 
A — > oo limit P However, we have now proved that its 
electrostatic energy (electron-background attraction plus 
background-background repulsion) is an approximation 
to the electrostatic potential of the exchange correlation 
hole, Eq. (12). This concept is further clarified in the 



Appendix, where the case of the uniform electron gas at 
low density is treated explicitly. 

The simplest approximation to the PC cell is a sphere 
of uniform density p(r) around the electron at position r 

with a radius r s (r) — (4^p(r)) fixed by the condition 
that the fictitious positive background exactly neutralizes 
the electron at its center. This leads to the simple PC- 
LDA approximation^ 



7/) LDA 



(r) = - 



10 



4tt\ 



1/3 



p(r 



,1/3 



(27) 



If we approximate the dipole moment of the cell in terms 
of the gradient of the density and we set it equal to zero 
we obtain the PC-GGA expression^ 



GGA 



(r) 



LDA 



(r) 



3 
350 



3 

-in 



1/3 



\Vp_(r) 
p(r 



1 7/3 



(28) 



In Fi g. [T] we compare the exact A —> 00 energy densities 
of Eq. J22T with the PC-LDA and PC-GGA approxima- 



tions of Eqs. ( 27 )-( 28 ) for the He atom, the spherical- 
ized B and C atoms and for the Ne atom, usi ng accur ate 
Hylleras and quantum Monte Carlo densities! 47 * 4 ^ 5 ! We 
see that the PC model becomes a rather good approxi- 
mation in the valence region of B, C and Ne, while being 
quite poor in the core region, and especially at the nu- 
cleus. The PC-LDA energy density is actually a better 
local approximation except close to the nucleus. The 
PC-GGA performs better globally (see Table [I]), but we 
clearly see that this is due to an error compensation be- 
tween the core region and the intershell region. 

The approximations done in the PC model are i) ne- 
glecting the cell-cell interaction, and ii) the gradient ex- 
pansion of Eqs. ( 27 )-( 28 1 which assumes a slowly varying 



density. At the nucleus, we can easily construct what 
would be the "exact" PC cell, so that we can at least 
remove approximation ii) and check the effect of approx- 
imation i) alone. The "exact" PC cell around the nucleus 
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FIG. 1: Energy density in the gauge of the electrostatic poten- 
tial of the exchange-correlation hole, Eq. ( 12 1, in the A — > 00 



limit. The exact SCE result of Eq. J22k is compared with the 



PC-LDA and PC-GGA approximations of Eqs. (27 1-1 28 1 
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TABLE I: Global value W^p] = f p(r)io 00 [p](r)dr for small 
atoms at different levels of approximation. The SCE cor- 
responds to the exact value, Eq. p2l, while PC-LDA and 



PC-GGA correspond, respectively, to Eqs. (271 and (28 1 
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C 


-4.7727 


-7.1446 


-5.6050 


Ne 


-8.0276 


-12.119 


-9.5463 



TABLE II: Comparison of the values at the nucleus of the 
energy density in the gauge of the exchange-correlation hole 
potential in the strong-interaction limit for small atoms. The 
value wsce(?~ = 0) corresponds to the exact expression of 

0) is the PC gradient expansion 



Eq. p2), the value Wpc A (r 



approximation of Eqs. <|27]) - ([28]) (the PC-LDA and PC-GGA 
are equal at the nucleusJTand wpc(r = 0) is the value from 
the "exact" PC cell of Eq. (30 1. 



is the sphere Qi of radius ai, with 



47rr p(r) dr = 1, 



and the "exact" value of wpc(r 
P(r) , 1 



0) is 



wpc(r = 0) = - 



dr- 



dr / dr 



(29) 



In Table [TT] we compare the values at the nucleus from 
the exact SCE, the PC-LDA or PC-GGA (they become 
equal at the nucleus) and the result of Eq. ( 30 ) for several 
atoms. We see that Eq. (30) is very accurate for N = 2 
electrons: in this case, in fact, when the reference elec- 
tron is at the nucleus, the other one is at infinity, so that 
the cell-cell interaction becomes indeed zero. For N > 2, 
we see that the "smearing hypothesis," i.e. the idea that 
the cell-cell interaction is negligeable, leads to some er- 
rors, although there is an improvement with respect to 
the gradient expansion of Eqs. ( 27 )-{ 28 1, reducing the 
relative error of about a factor 2.^1ong these lines, one 
might try to construct an improved PC model that per- 
forms locally better than the PC-GGA, which, as said, 
achieves good global accuracy at the price of error com- 
pensation between different regions of space. 



V. ENERGY DENSITIES ALONG THE 
ADIABATIC CONNECTION 

A. Kohn-Sham (A = 0) 

At zero coupling-strength the exact solution for the 
wave function becomes a Slater determinant $ = 
|0i...0at), and the energy density iuo[p]( r ) in the gauge of 
the exchange-correlation hole is given by the electrostatic 
potential of the KS exchange hole h x (r, r'), 



If h x (r,r') , 
-oM(r) = - /— -dr, 



as the pair density simply writes 

P 2 °(r,r') =p(r)p(r')+p(r)h x (r,r'). 



(31) 



(32) 



One can use in Eq. ( |31[ ) the exact exchange hole built 
from a Hartee-Fock like expression in terms of the KS 
orbitals </>i, or a density functional approximation for 
h x (r,r'), e.g., the one of Becke and RousselP^I These 
two choices would correspond, respectively, to construct 
a hyper-GGA and a meta-GGA functional from a local 
interpolation along the adiabatic connection. 

The aim of the present work is a preliminary study 
of exact energy densities along the adiabatic connec- 
tion. The exact KS orbitals and the corresponding non- 
interacting potential V® xt for a given physical density can 
be found in an exact way, e.g,. by inversion of the KS 
equationj^H^o r by the use of Lieb's Legendre transform 
DFT formalismPMMl 

For an ISI-like interpolation on the energy density, 
u>o[p](r) will be a key ingredient. Additionally, knowledge 
of the next leading order in the asymptotic expansion of 
the local energy density around A = is necessary, but 
not available yet. The next leading order in the asymp- 
totic expansion constitutes an active field of research in 
our group (see also the discussion in Sec. VII). 



B. Physical (A = 1) 

To compute the exact energy density at coupling- 
strength A = 1 we resort to 



wi[p](r) = 



1 



2p(r) 



|r-r'| 2 



P(r') 
|r-r' 



dr', 
(33) 

with the pair density given by the full many-body wave 
function in Eq. (14). The density pi(r) out of 
P 2 1 (r,r / ) defines the density p(r) = pi(r) to be held con- 
stant along the adiabatic connection. 

The exact wi[p](r) can serve as benchmark for models 
on w\ [p] (r) but also gives an estimate of the importance 
of the strong- interaction limit in the iux[jo](r) model. If 
the physical system is close to the KS one, correlation 
is less important and already Hartree-Fock should per- 
form well. In this case, we expect that inclusion of the 
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A — > oo information in the w\ [p] (r) model does not lead 
to a major improvement. In contrast, for stronger cor- 
related system the physical energy density should tend 
more towards the A — > oo limit and the SCE concept can 
provide useful input for an accurate model for w>[p](r). 
The relevance of the strong-interaction limit will be dis- 
cussed in the next section. 



C. Results 

1. Coulomb external potential 

We have performed full-CI calculations in an aug-cc- 
pVTZ basis for some two and four electron atoms within 
the Gamess-US package^ to obtain an accurate ground 
state wave function for the physical interaction strength. 
Starting from this, we are able to calculate the energy 
density in the gauge of the exchange-correlation hole for 
A = 0, 1, oo. 

At A = 1 we calculate the energy density from the full- 
CI pair density, Eq. (33 ), by a program similar to the one 

inRef.Hl 

we have to 



(31) 



used for the calculation of w con d 

For the energy density at A = 0, Eq. 
compute the single particle KS orbitals corresponding to 
the full-CI density first. In the case of 2 electron atoms 
they are readily constructed by the simple relation 



0(r) 



P(r) 
2 



(34) 



For the four electron atoms we cho ose th e scheme of 
van Leeuwen, Baerends and Gritsenkd 51 * 52 ! to invert the 
KS equations. In the strong-interaction limit we calculate 
the energy density within the SCE concept, see Sec. |IV A"| 
and Refs. US 

In Fig. [2] we show the energy densities for A = 0, 1, oo 
for two- and four-electron atoms. As expected, He and Be 
are relatively weakly correlated, and their A — 1 energy 
densities are much closer to the KS ones than to the 
SCE ones. Here a description at the Hartree-Fock level 
is very reasonable and gives indeed at least 98.5% of the 
total energy. The anion H~, instead, being a system with 
a more diffuse density and thus more correlated, has a 
physical energy density that is much more in between the 
KS and the SCE curves, with an Hartree-Fock treatment 
giving only 94% of the total energy. Here we expect the 
inclusion of the information from the strong-interaction 
limit to be important. The valence regions of Be and Li~ 
(see the insets in Fig. [2| can also be better described by 
a proper inclusion of the A = oo information. 



~ -0.8 




0.5 1 1.5 2 2.5 3 3.5 
distance from nuclei (a Bohr ) 






-0.5 
-1 

-1.5 
-2 

-2.5 



1 2 3 4 5 6 7 
distance from nuclei (a Bohr ) 




12 3 4 

distance from nuclei (a Bohr ) 



-0.07 
-0.14 E 
-0.21 
-0.28 



i 4 
X = 
X=\ 

X — ¥ oo 



1 2 3 4 5 
distance from nuclei (a Bohr ) 



2. Harmonic external potential 

Another useful class of systems to investigate the im- 
pact of the strong-interaction limit on the physical en- 
ergy density is given by model quantum dots, where elec- 



FIG. 2: Energy densities in the gauge of the electrostatic 
potential of the exchange-correlation hole for accurate full-CI 
densities (aug-cc-pVTZ) and coupling strength A = 0, 1 and 
oo. 
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FIG. 3: Energy densities in the gauge of the electrostatic 
potential of the exchange-correlation hole for a quantum dot 
with N — 2 electrons with less pronounced correlation (u = 
0.5) and pronounced correlation (uu = 0.0014). 



trons are confined by a harmonic potential and correla- 
tion gains importance as the spring constant is lowered. 
We have computed the energy density for N = 2 elec- 
trons in three dimensions for spring constants for which 
an analytic solution for the wave function can be foundl^l 
The results are displayed in Fig. [3] for the largest and 
smallest spring constant considered. As expected, the 
physical energy density comes closer to the SCE energy 
density in the stronger correlated case. Additionally, a 
remarkable feature we observe is that the physical energy 
density crosses the SCE energy density. By intuition one 
would expect the physical energy density to be always in 
between the KS and SCE energy densities, as the KS en- 
ergy density represents the weakest possible correlation 
and the SCE energy density the strongest possible corre- 
lation in the given density. However, the wave functions 
are chosen according to the global quantities 



min(*|T|*) 

-t>—yp 

mm(^\V ee \^) 



min(\f/|T 



v ee m 



4-KS 

^SCE 

*A=1, 



(35) 
(36) 
(37) 



FIG. 4: Energy density in the gauge of the electrostatic po- 
tential of the exchange-correlation hole for a N = 2-electron 
model quantum dot calculated with the full KS, physical and 
SCE pair density and with the asymptotic expansion of the 
physical pair density, Eq. | |39[ ). The KS and 1st order curve 
lie on top of each other. 



yielding the global inequalities 



(*SCE|K e |*SCE> < (*A=l|Ke|*A=l) < (*KS \V ee | * KS ) ■ 

.(38) 

Locally, these inequalities can be violated without violat- 
ing the global ones, and hence the physical energy density 
can go below the SCE energy density. 

The crossing feature can be attributed to polarization 
effects, that are present in the physical case, but are not 
contained in the KS and SCE pair density. In the KS case 
we rely on the independent-particle picture and polariza- 
tion of a particle due to other particles is not described 
by this model. Whereas for the SCE, we build on point 
particles that are a priori unpolarizable. To underline 
this argument we have computed the asymptotic behav- 
ior of the physical energy density for the quantum dot 
with ui = 0.5 by use of the asymptotic expansion of the 
physical pair density^ 



P(r,r') 
p(r)pjv-i(r') 



1 



•2 — cos(Aft) + - 
r 



(r -> oo), (39) 



where /9jv_i(r') is the density of the N— 1-particle system 
and Af2 the angle between r and r'. The second order 
term in Eq. ( 39 ) represents the polarization correction. 
As can be seen from Fig.[4j for large distances |r — r'| the 
energy density in the KS and SCE case behaves like ~ 1 /r, 
which corresponds to the first order in the expansion of 
the physical pair density. Inclusion of the second order 
term in the energy density gives essentially the physi- 
cal behavior and deviates from the KS and SCE energy 
densities. 

Although the crossing happens in a region in which 
the density is very small and thus with an almost negli- 
gible energetic contribution, the analysis presented here 
can be helpful in constructing models for w\[p](r). No- 
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tice that, instead, with the Coulomb external poten- 
tial we always observed, so far, the expected behavior 
WA^oo(r) < w\=i(r) < w x =o(r) everywhere. 



THE LOCAL FORM OF THE 
LIEB-OXFORD BOUND 



VI. 



The Lieb-Oxford (LO) bouncP^fi] j s a rigorous lower 
bound to the indirect part of the electron-electron repul- 
sion energy VF[ , I / ] associated with a given many-electron 
wave function "J, 



W[V] = (VlVee]*) - Ufa] > -C / drp*(r) 



n4/3 



(40) 



where /c*(r) is the density obtained from the wave func- 
tion ^. Th e posi tive constant C is rigorously kno wn to 
have a valueP122l q < 1.679. It has been suggestecP^ 
that a tighter bound can be obtained by taking the value 
of C that corresponds to the low-density limit of the uni- 
form electron gas, C « 1.44, since the bound is known 
to be more challenged when the number of electrons 
increases^ and when the system has low density.^ 

The LO bound translates into a lower bou nd for the 
exchange and exchange-correlation functionals j 61 l 63 l 



E x [p]>E xc [p}>-C J drp(r) 4 / 3 , 



(41) 



simply because E. x [p] = Wa=o [p] is the indirect Coulomb 
repulsion of the Slater determinant of KS orbitals, and 
E xc [p] is the sum of the indirect Coulomb repulsion of 
the physical wave function, Wa=i[/>], plus the correlation 
correction to the kinetic energy which is always positive. 

The way the LO bound is used in the construction of 
approximate functionals is, usually (with the exception 
of Ref. [Ml), by imposing it locally (see, e.g., Refs. \TT[ 
I23p . That is, a given approximate exchange-correlation 
functional, EWf[p] = J p(r) e^Ff (r) dr, is required to 
satisfy 



(42) 



This is a sufficient condition to ensure the global bound 
of Eq. (41 1, but it is by no means necessary (see, e.g., 
Ref. In other words, there is no proof that a lo- 

cal version of the LO bound should hold. Actually, be- 
fore even asking wether a local version of the LO bound 
should hold or not, we need to understand to which def- 
inition (gauge) of the energy density should apply the 
local LO bound of Eq. (42). In fact, since energy densi- 



ties are not uniquely defined, the inequality (42 1 should 



be satisfied only for a well defined gauge: one can indeed 
always add to e^^(r) a quantity that integrates to zero 
and violates Eq. ( |42[ ) in some region of space. 

We argue here that i) the gauge of the local LO bound 
is the conventional one of the electrostatic energy of 
the exchange-correlation hole, and ii) that the local LO 



bound is certainly violated, at least in the tail region of 
an atom or of a molecule, and in the bonding region of 
a stretched molecule. The argument behind point i) is 
the following. For a given density p, the wave function 
^[p] that maximally challenges^ the LO bound is the 
one that minimizes the expectation ($[/o]|V ee | 'I' [/>]), i.e., 
by definition, 'J'sce [/£>]• In fact, we also have 



E x [p] > E xc [p] > Wnlp] > -C j dvp{i 



,4/3 



(43) 



In Sec. |IV| we have discussed the energy density associ- 
ated to Woo [p] in the gauge of the electrostatic potential 
of the exchange-correlation hole. We have also shown 
that this energy density can be approximated by the PC 
model that considers the electrostatic energy of a cell 
around the reference electron of positive charge p+(r) = 
p(r) . The LDA version of this approximation has exactly 
the same form of the local LO bound. Moreover, the 
recently suggested valued C « 1.44 is extremely close 
to the one of the PC-LDA model, C PC ~ 1.45. Notice 
that the fact that the PC model is in the gauge of the 
electrostatic energy of the exchange-correlation hole fol- 
lows from the properties of the strong-interaction limit 
of DFT, in particular Eq. (l25|). If the PC model is an ap- 



proximation in this gauge, and if the LO bound is locally 
equal to it, then conclusion i) should follow. 



We then easily see that the local LO bound of Eq. ( 42 ) 



is certainly violated in the tail region of an atom or a 
molecule, where the exact energy density in the con- 
ventional exchange-correlation hole gauge goes like — 1/r 



while the right-hand side of Eq. ( 42 1 decays exponen- 
tially. The local bound is also violated in the bond re- 
gion of a stretched molecule. As an example, we show in 
Fig. [5] the energy densities of the stretched H2 molecule 
for A = and A = 1: with C = 1.44 the local bound is vi- 
olated in the bond region when the intcrnuclcar distance 
is R > 7 a.u., and with C = 1.67 when R>8 a.u. 

As a concluding remark, we can say that it is very 
difficult, or maybe even impossible, to find a rigorous 
local lower bound for the energy density. In fact, we 
have just seen in Sec. [Vjthat, at least for the harmonic 
external potential, it is not even true that WA=i(r) > 
Woo(r) everywhere. This means that even if we maximize 
the correlation between the electrons we do not construct 
a rigorous local lower bound, but only a global one. 



VII. CONCLUSIONS AND PERSPECTIVES 

We have derived an exact expression for the energy 
density in the strong-interaction limit of DFT in the 
gauge of the exchange-correlation-hole electrostatic po- 
tential, and we have computed it for small atoms and 
model quantum dots. A careful analysis of the point- 
charge plus continuum (PC) model showed that this ap- 
proximation is formulated in the same gauge, and a com- 
parison with the exact results showed that it is locally 
reasonable. 
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a local slope at A = is to produce exact results for 
this quantity, crucial to assess approximations. This can 
be achieved with the Legendre transform techniques de- 
veloped in Refs. I29I30I and it is the object of a current 
project. A possible way, then, to construct an approxi- 
mate local slope is to use the so-called "extended Over- 
hauser model"'^HS2l locally, in a perturbative way. A local 
next leading term at A — oo can also be constructed by 
deriving the exact exchange-correlation hole correspond- 
ing to the wave function of the zero-point oscillations, 
discussed in Ref. [9l All this will be investigated in future 
work. 
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FIG. 5: Violation of the local form of the Lieb-Oxford bound 
for the stretched H2 molecule. 



Our formalism also allowed to analyze the local ver- 
sion of the Lieb-Oxford bound: we were able to assess to 
which gauge should the local LO bound correspond, then 
showing that it is certainly violated. Our findings are in 
agreement with (and give formal support to) the very 
recent results of Vilhena et alW^ (which only appeared 
when this manuscript was completed). More generally, 
our results suggest that it is very difficult (if not impossi- 
ble) to derive a rigorous local lower bound for the energy 
density. 

We have also discussed the idea of a local interpolation 
along the adiabatic connection. The values of the local 
energy density in the same gauge at A = and A = 00 
are now available, cither exactly or in an approximate 
way. Even if we have found that in the harmonic exter- 
nal potential the physical energy density is not always 
in between the A = and the A = 00 curves, the re- 
gions of space in which the expected order is reversed are 
energetically not important. In the external Coulomb 
potential we have found, instead, the expected behavior 
WA^oo(r) < w\=x(r) < wa=o(i") everywhere. 

To really be able to build a local interpolation, at least 
the slope at A = 0, and possibly the next leading term 
at A = 00, are also needed in a local form, and in the 
same gauge. A first step towards the construction of 



Appendix A: PC cell and xc-hole 

In this Appendix we clarify the difference between the 
exchange-correlation hole and the PC cell by consider- 
ing the uniform electron gas in the extreme low-density 
limit, further extending the argument already given in 
the Appendix of Ref. |H 

More than seventy years ago, Wignei^HU pointed out 
that electrons embedded in a compensating uniformly 
charged background would crystallize at sufficiently low 
values of the density p. The SCE construction can be 
seen as nothing else than the Wigner idea generalized to 
a nonuniform density p(r). Indeed, in Ref. M the SCE 
formalism is presented as a "floating" Wigner crystal in 
a non-euclidean space, with the metric determined by the 
density p(r). 

In the case of the uniform electron gas, the SCE co- 
motion functions are simply the positions of the bec 
lattice points with origin fixed at the reference elec- 
tron. Notice that the constraint that the density is 
uniform, Eq. (36), forces us to consider a "floating" 



Wigner crystal, which corresponds to the linear super- 
position of all the possible origins and orientations of the 
crystal, thus restoring the translational symmetry. The 
exchange-correlation hole p(g(r) — 1), with g(r) the pair- 
distribution function, can then be simply constructed by 
considering that the expected numbers of electrons in a 
spherical shell of radius r around the reference electron 
at the origin is given by 

dN(r\0) = pg(r)4:7rr 2 . (Al) 
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(perfect long-range order). Notice that the exchange- 
correlation hole for the broken symmetry solution (with- 
out translational invariance) would be, instead, much less 
structured, but here we are interested in the solution con- 
strained to the uniform density. The way the electrostatic 
energy is calculated from the PC cell and the exchange- 
correlation hole is also different P 
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FIG. 6: PC cell divided by the density p (i.e., c(r) — —9(r„ — 
r)), and the pair-correlation function g(r) — 1, corresponding 
to the exchange-correlation hole divided by the density p for 
the extreme low-density electron gas. 



We can then place very narrow normalized gaussians (al- 
most delta functions) at the bec sites around the refer- 
ence electron and take the spherical average. This way, 
we obtain the extreme low-density limit of g{r). In Fig. [6] 
we compare this low-density (or SCE) g(r) — 1 with the 
PC cell c(r) in the same units, c(r) = —9{r s — r), with 
0(x) the Heaviside step function. We see that the two 
are very different, except for r/r s < 1. The exchange- 
correlation hole has positive peaks (indicating the po- 
sitions of the other electrons) that extend to r — > oo 



When we use t he e xchange-correlation hole to evaluate 
the energy, Eq. ( A2 ) , we need to evaluate an infinite sum 
(all the peaks in Fig. [6]) which converges very badly (the 
Madelung sum), and that can be dealt with, for example, 
the Ewald method. When we use the PC cell, instead, 
we face two very simple, short-ranged, integrals.^ The 
two results differ only by 0.45%, as it was already noted 
in Ref. [72j where it was also proven that the PC value 
is a rigorous lower bound for the energy of the uniform 
electron gas. 

Notice that if, instead, we consider the PC cell as 
a model for the exchange correlation hole and we use 
c(r) in Eq. (A2| instead of g(r) — 1, we get a very poor 
resullP with an error of ~ 17%. The PC model does 
approximate the electrostatic potential of the exchange- 
correlation hole by constructing it in a different way. 
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